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Abstract. In this paper we consider disjoint decomposition of algebraic and non-linear par- 
tial differential systems of equations and inequations into so-called simple subsystems. We 
exploit Thomas decomposition ideas and develop them into a new algorithm. For algebraic 
systems simplicity means triangularity, squarefreeness and non-vanishing initials. For dif- 
ferential systems the algorithm provides not only algebraic simplicity but also involutivity. 
The algorithm has been implemented in Maple. 



1 Introduction 

Nowadays, triangular decomposition algorithms, which go back to the characteristic set method 
by Ritt [Rit50] and Wu [WuOO], and software implementing them have become powerful tools for 
investigating and solving systems of multivariate polynomial equations. In many cases these meth- 
ods are computationally more efficient than those based on construction of GrObner bases. As an 
example of such problems one can indicate BoOLean polynomial systems arising in cryptanalysis 
of stream ciphers. For those systems triangular decomposition algorithms based on the character- 
istic set method revealed their superiority over the best modern algorithms for the construction 
of Grobner bases [sGH09]. 

For terminology, literature, definitions and basic proofs on triangular-decomposition algo- 
rithms for polynomial and differential-polynomial systems we refer to the excellent tutorial papers 
[Hub03a, Hub03b] and to the bibliographical references therein. 

Among numerous triangular decompositions the Thomas one stands by itself. It was suggested 
by the American mathematician J.M.Thomas in his books [Tho37, Tho62] and decomposes a finite 
system of polynomial equations and inequations into finitely many triangular subsystems that he 
called simple. Unlike other decomposition algorithms it yields a disjoint zero decomposition, that 
is, it decomposes the affine variety or quasi-affine variety defined by the input into a finite number 
of disjoint quasi-affine varieties determined by the output simple systems. Every simple system is 
a regular chain. 

During his research on triangular decomposition, Thomas was motivated by the Riquier- 
Janet theory [RiqlO, Jan29], extending it to non-linear systems of partial differential equations. 
For this purpose he developed a theory of (Thomas) monomials, which generate the involutive 
monomial division called THOMAS division in [GB98a]. He gave a recipe for decomposing a non- 
linear differential system into algebraically simple and passive subsystems [Tho37] . 

Differential Thomas decomposition differs noticeably from that computed by the famous 
Rosenfeld-Grobner algorithm [BLOP09, BLOP95] which forms a basis of the diffalg and BLAD 
libraries [BH04, Bou09] as well as from other differential decompositions (e.g. [BKRM01]). We 
found that diffalg and BLAD are optimized and well-suited for ordinary differential equations. 
However, all other known methods give a zero decomposition which, unlike that in THOMAS de- 
composition, is not necessarily disjoint. 

A first implementation of the THOMAS decomposition was done by Teresa Gomez-Diaz in AX- 
IOM under the name "dynamic constructible closure" which later turned out to be the same as the 
Thomas decomposition [DelOO]. Wang later designed and implemented an algorithm construct- 
ing the Thomas decomposition [Wan98, WanOl, LW99]. For polynomial and ordinary differential 
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systems Wang's algorithm was implemented by himself in Maple [Wan04] as part of the software 
package epsilon [Wan03], which also contains implementations of a number of other triangular 
decomposition algorithms. A modified algorithmic version of the Thomas decomposition was con- 
sidered in [Ger08] with its link to the theory of involutive bases [GB98a, Ger05, Ger99]. The latter 
theory together with some extensions is presented in detail in the recent book [SeilO]. 

In the given paper we present a new algorithmic version of the Thomas decomposition for 
polynomial and (partial) differential systems. In the differential case the output subsystems are 
Janet involutive in accordance to the involutivity criterion from [Ger08], and hence they are 
coherent. Moreover, for every output subsystem the set of its equations is a minimal Janet basis 
of the radical differential ideal generated by this set. The algorithm has been implemented in 
Maple for both the algebraic and differential case. For a linear differential system it constructs a 
Janet basis of the corresponding differential ideal and for this case works similarly to the Maple 
package Janet (cf. [BCG+03]). 

This paper is organized as follows. In §2 we sketch the algebraic part of our algorithm for the 
Thomas decomposition with its main objects defined in §2.1. The algorithm itself together with 
its subalgorithms is considered in §2.2. Decomposition of differential systems is described in §3. 
Here we briefly introduce some basic notions and concepts from differential algebra (§3.1) and 
from the theory of involutive bases specific to Janet division (§3.2) together with one of the two 
extra subalgorithms that extend the algebraic decomposition to the differential one. The second 
such subalgorithm is considered in §3.3 along with the definition of differential simple systems. 
Subsection §3.4 contains a description of the differential Thomas decomposition algorithm. Some 
implementation issues are discussed in §4, where we also demonstrate the Maple implementation 
for the differential decomposition using the example of a system related to control theory. 

We omit the proofs for compactness. They will be published elsewhere. 

2 Algebraic Systems 

The algebraic Thomas decomposition deals with systems of polynomial equations and inequations. 
This section introduces the concepts of simple systems and disjoint decompositions based on 
properties of the set of solutions of a system. A pseudo reduction procedure and several splitting 
algorithms on the basis of polynomial remainder sequences are introduced as tools for the main 
algorithm, which is presented at the end of the section. 

2.1 Preliminaries 

Let F be a computable field of characteristic and R := F [x±, . . . , x n ] the polynomial ring 
in n variables. A total order < on the indeterminates of R is called a ranking. The notation 
R = F[xi, . . . , x n ] shall implicitly define the ranking x\ < . . . < x n . The indeterminate x is called 
leader of p G R if x is the <-largest variable occurring in p and we write ld(p) = x. If p G F, we 
define ld(p) = 1 and 1 < x for all indeterminates x. The degree of p in ld(j>) is called rank of p 
and the leading coefficient init(p) G F[ y \ y < ld(p) ] of \d(p) lank ^ in p is called initial of p. 

For a G F , where F denotes the algebraic closure of F, define the following evaluation 
homomorphisms: 



For a polynomial p G R, the symbols p = and p^ shall denote the equation p = and inequation 
p ^ 0, respectively. A finite set of equations and inequations is called an (algebraic) system over 
R. Abusing notation, we sometimes treat p = or p-t as the underlying polynomial p. A solution of 
a system 5 is a tuple a G F such that <ft a (p) = for all equations p— G S and <^ a (p) ^ for all 
inequations p^ G S. The set of all solutions of S is denoted by &oi{S). 



a : F[xi , . . . , x n ] -> F : x,- L i-> a, t 
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Define S x := {p € S | ld(p) = x}. In a situation where it is clear that |5 X | = 1, we also use 
S x to denote the unique element of S x . The subset S <x := {p G S | ld(p) < x} can be considered 
a system over F[ y \ y < x ]. Furthermore, the sets of all equations p— G S and all inequations 
p^ G S are denoted by S = and , respectively. 

The general idea of the THOMAS methods is to use the homomorphism 4>< x , a to treat each poly- 
nomial p G S x as the univariate polynomial 4>< x , a {p) G F[x] for all a G &ol(S<: x ) simultaneously. 
This idea forms the basis of our central object, the simple system: 

Definition 2.1 (Simple Systems). Let S be a system. 

1. 5" is triangular if \S Xi \ < 1 V 1 < i < n and Sfl {c=,cy | c G F} = 0. 

2. £ has non-vanishing initials if (/> a (init(p)) ^ V a G 6o[(S< Xi ) and p G S Xi for 1 < i < n. 

3. S* is square-free if the univariate polynomial <p< Xi ,a{p) G F[xi] is square-free V a G &ol(S <Xi ) 
and p G S Xi for 1 < i < n. 

4. 5* is called simple if it is triangular, has non-vanishing initials and is square-free. 

Although all required properties are characterized via solutions of lower-ranking equations and 
inequations, the Thomas decomposition algorithm does not calculate solutions of polynomials. 
Instead, it uses polynomial equations and inequations to partition the set of solutions of the 
lower-ranking system to ensure the above properties. 

Remark 2.2. Simplicity of a system guarantees the existence of solutions: If b G &ol(S <x ) and 
S x is not empty, then 4>< x .h(S x ) is a univariate polynomial with exactly r&nk(S x ) distinct roots. 
When extending b to a solution (b,a) of S< Xl for an equation in S x there are rank(S , a; ) choices 
for a, whereas for an inequation or empty S x all but finitely many a G F give an extension. 

To transform a system into a simple system, it is in general necessary to partition the set 
of solutions. Instead of an equivalent simple system, this leads to a so-called decomposition into 
simple systems. 

Definition 2.3. A family (5 , ,)™ 1 is called decomposition of S if 6oI(5) = [JiLi 6o K S i)- A 
decomposition is called disjoint if &ol(Si) n &ol(Sj) = V i ^ j. A disjoint decomposition of a 
system into simple systems is called (algebraic) THOMAS decomposition. 

For any algebraic system S, there exists a Thomas decomposition (cf. [Tho37], [Tho62], 
[Wan98]). The algorithm presented in the following section provides another proof of this fact. 
First, we give an easy example of a THOMAS decomposition. 

Example 2.4. Consider the equation 

p = y — x r — x . 
A Thomas decomposition of {p=\ is given by: 

({(y 2 x 3 x 2 ) = 1 (x-(x + 1))*}, {y=, (x ■ (x + 1))=}) 

2.2 Decomposition Algorithms 

Our version of the decomposition algorithm in each round treats one system, potentially splitting 
it into several subsystems. For this purpose, one polynomial is chosen from a list of polynomials 
to be processed. This polynomial is pseudo-reduced modulo the system and afterwards combined 
with the polynomial in the system having the same leader. To ensure that all polynomials are 
square-free and their initials do not vanish, the system may be split into several ones by initials 
of polynomials or subresultants. 

From now on, a system 5* is presented as a pair of sets (St, Sq), where St represents a candidate 
for a simple system while Sq is the queue of elements to be processed. St is always triangular and 
(St)x denotes the unique equation or inequation of leader x in St, if any. St also fulfills a weaker 
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form of the other two simplicity conditions, i.e., for any solution a of (St)< x U {Sq)<x, we have 
</> a (init((S , T)x)) and <j> <X!a {{S T )x) is square-free. 

From now on, let prem be a pseudo remainder algorithm 3 in R and pquo the corresponding 
pseudo quotient algorithm, i.e., for p and q with ld(p) = ld{q) = x 

m ■ p = pquo(p, q, x) ■ q + prem(p, q, x) (1) 

where deg x (g) > deg a ,(prem(p, q, x)) and m £ R\ {0} with ld(m) < x and m | init(q , ) fc for some 
k £ Z>o- Note that if the initials of p and q are non-zero, the initial of pquo(p, q, x) is also non-zero. 
Equation (1) only allows us to replace p with prem(p, q,x) if m does not vanish on any solution. 
The below Algorithm (2.5) and Remark (2.6) require the last property, which, by definition, holds 
in simple systems. 

The following algorithm employs pseudo remainders and the triangular structure to reduce a 
polynomial modulo St' 

Algorithm 2.5 (Reduce). 

Input: A system S, a polynomial p £ R 

Output: A polynomial q with 4> a (p) = if and only if (f) a (q) — for each a £ &ol(S). 
Algorithm: 
1: X 4— ld(p); q 4- p 

2: while x > 1 and (St) x is an equation and rank(g) > rank ((St) a;) do 
3: q 4- prem((7, (S T ) x ,x) 
4: x <- \d(q) 
5: end while 

6: if x > 1 and Reduce(5, init(g)) = then 
7: return Reduce^, q - init(g)x rank ^)) 
8: else 

9: return q 
10: end if 

A polynomial p is called reduced modulo St if Reduce^, p) = p. A polynomial p reduces to q 
modulo St if Reduce(S',p) = q. 

The result of the Reduce algorithm does not need to be a canonical normal form. It only needs 
to detect polynomials that vanish on all solutions of a system: 

Remark 2.6. Let p £ R with ld(p) = x. Reduce(5,p) = implies (f> a (p) = V a £ &ol(S< x ). 

The converse of this remark only holds if (Sq)< x = 0, i.e., (St)< x is simple. If it is not simple, 
but ld(p) = x and (Sq)^ x — hold, we still have some information. In particular, Reduce(5, p) ^ 
implies that either 6o[(6' <x ) is empty or there exists a £ &o\(S <x U {{St) x }) such that </> a (p) 7^ 0. 

We now direct our attention to the methods we use to produce disjoint decompositions. Since 
(S U {p^} , S U {p=}) is a disjoint decomposition of S, we will use the following one-line subalgo- 
rithm as the basis of all the splitting algorithms described below. 

Algorithm 2.7 (Split). Input: A system S, a polynomial p £ R 
Output: The disjoint decomposition (S U {p^} , S U {p=}) of S. 
Algorithm: 

1: return {{St, Sq U {p^}) , {St, Sq U {p=})) 

The output of the following splitting algorithms is not yet a disjoint decomposition of the input. 
However, the main algorithm Decompose will use this output to construct a disjoint decomposition. 
We single out these algorithms to make the main algorithm more compact and readable. For details 
we refer to the input and output specifications of the algorithms in question. 

The algorithm InitSplit ensures that in one of the returned systems the property 2 in Definition 
(2.1) holds for the input polynomial. In the other system the initial of that polynomial vanishes. 

3 In our context prem does not necessarily have to be the classical pseudo remainder, but any sparse 
pseudo remainder with property (1) will suffice. 
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Algorithm 2.8 (InitSplit). Input: A system S, an equation or inequation q with \d(q) = x. 
Output: Two systems Si and S 2 , where (Si U {q},S 2 ) is a disjoint decomposition of 5* U {q}. 
Moreover, </> a (init(g)) ^ holds for all a 6 60 [(Si) and </> a (init(g)) = for all a G &ol(S 2 )- 
Algorithm: 

1: {S U S 2 ) ^Split(5,init(g)) 

2: if q is an equation then 

3: {S 2 )q <" (5 2 )Q U {(« - init(g)x rank (")) = } 
4: else if q is an inequation then 

5: {S 2 )q (&)q U {( 9 - init(g)x- nk (^)^} 
6: end if 

7: return {Si,S 2 ) 

In Definition (2.1) we view a multivariate polynomial p as the univariate polynomial </><id( P ),a(p)- 
For ensuring triangularity and square- freeness, we often compute the gcd of two polynomials, which 
generally depends on the inserted value a. Subresultants provide a generalization of the ElJCLlDean 
algorithm useful in our context and their initials distinguish the cases of different degrees of gcds. 



Definition 2.9. Let p,q G R with ld(p) = ld(g) = x, deg x (p) = d p > deg x (q) = d q . We denote by 
PRS(p, q, x) the subresultant polynomial remainder sequence (see [Hab48], [Mis93, Chap. 7], 
[YapOO, Chap. 3]) of p and q w.r.t. x, and by PRSi(p, q, x), i < d q the regular polynomial of degree 
i in PRS(p, q,x) if it exists, or otherwise. Furthermore, PRSd p [p, q, x) := p, PRSd q (p,q,x) := q 
and PRSi(p, q, x) := 0, d q < i < d p . 

Define resj(p, q,x) := init (PRSj {p, q, x)) for < i < d p , whereas resd p (p, q, x) := 1 and 
reso(p, q, x) := PRSo (p, q, x). Note that res(p, q, x) := resrj(p, q, x) is the usual resultant. 

Definition 2.10. Let S be a system and pi,p 2 G R with ld(pi) = ld{p 2 ) = x. If |©o[(S< x )| > 0, 
we call 



i := min{z G Z> | 3 a e &ol(S <x ) such that deg a ,(gcd(^ <a; , a (pi), 4> <x . a (p 2 ))) = i} 

the fiber cardinality of p-y and p 2 w.r.t. S. Moreover, if {Sq)< x = 0, then 

i' := min{z G Z>o | Reduce(reSj(pi,p2, x), St) = V j < i and Reduce(resj(pi,p2, x), St) 7^ 0} 

is the quasi fiber cardinality of p\ and p 2 w.r.t. S. A disjoint decomposition {Si, S 2 ) of S such 
that 

1. deg a .(gcd(0 <Xia (pi),^ <Xia (p 2 ))) = iVae &ol{{Si) <x ) 

2. deg x (gcd(^ <IiB (pi),0 <XiB (p2))) > i V a g Sol ((5 2 ) <a! ) 

is called the i-th fibration split of pi and p 2 w.r.t. S 1 . A polynomial re J? with ld(r) = x such 
that deg x (r) = i and 

4><xA T ) ~ gcd(0 <a; ,a(Pl), 0<x,a(P2)) V a G 6o[ ((Sl)<x) 

is called the «-th conditional greatest common divisor of pi and p 2 w.r.t. S, where p ~ q if 
and only if p £ K q. Furthermore, q G R with ld(q) = a; and deg x {q) = deg x {pi) — i such that 

4>< X AQ) ~ 1(1 ^Ti^J ( \\ V a G 6o( {(Si) <x ) 

gcd(0 <X;a (pi),0 <a;;a (p2)) 

is called the i-th conditional quotient of pi by p 2 w.r.t. S. By replacing a (p 2 ) in the above 
definition with J^(^< x ,a(f>i))i we get an i-th square-free split and i-th conditional square-free 
part of pi w.r.t. S. 
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The fiber cardinality is often not immediately available, as we may be unable to take inequations 
into account. However, we can use all information contained in the equations using reduction, if 
all equations are contained in St- Thus we require (Sq)< x = before doing any reduction. 

In this situation, the quasi fiber cardinality is easy to calculate and in many cases will be 
identical to the fiber cardinality. Furthermore, if we consider the system S2 from an i-th fibration 
split of some polynomials for a system S and ensure that ((S 2 )q)<,,. = ®, then the quasi fiber 
cardinality of the same polynomials for £2 will be i + 1. Therefore and due to the following lemma, 
the quasi fiber cardinality is good enough for our purposes. 

Lemma 2.11. Let \&ol(S <x )\ > and (Sq)< x = 0. For p\, p 2 as in Definition (2.10) with 
a (init(pi)) 7^ V a e &ol(S <x ) and rank(pi) > rank(p 2 ), let i be the fiber cardinality of p\ and 
pi w.r.t. S and i' the corresponding quasi fiber cardinality. Then 



where the equality holds if and only if \&ol(S <x U {res,/(pi,p2, z)^})| > 0. 

Corollary 2.12. Let |<5o[(S< x )| > and (Sq)< x = 0- For polynomials p\, p 2 as in Definition 
(2.10) with </> a (init(pi)) 7^ and </> a (init(p 2 )) 7^ V a 6 &ol(S <x ), let i be the fiber cardinality of 
Pi and pi w.r.t. S and i' the quasi fiber cardinality of pi and prem(p2,pi,x) w.r.t. S. Then 



with equality if and only if \&ol(S <x U {resj/(pi, prem(p 2 ,pi, x), x)^})\ > 0. 

The following algorithm calculates the quasi fiber cardinality of two polynomials. It is used as 
the basis for all algorithms that calculate a greatest common divisor or a least common multiple. 

Algorithm 2.13 (ResSplit). Input: A system S with (Sq)< x = 0, two polynomials p,q 6 R with 

\A{p) = ld(g) = x, rank(p) > rank(q) and </> a (init(p)) =/= for all a G &ol(S <x ). 

Output: The quasi fiber cardinality i of p and q w.r.t. S and an i-th fibration split (Si,S 2 ) of p 

and q w.r.t. S. 

Algorithm: 

1: i <— min{i e Z>o | Reduce(reSj (p, q 7 x), St) = V j < i and Reduce(resi(p, q, x), St) 7^ 0} 
2: return (i, Si, S2) := (i, Split(S', resi(p, q, x))) 

Similarly to the InitSplit algorithm (2.8), the following algorithm does not return a disjoint 
decomposition, but Decompose uses it to construct one. 

Algorithm 2.14 (ResSplitGCD). Input: A system S with (Sq)< x = 0, where (St)x is an equation, 
and an equation q = with \d(q) — x. Furthermore rank(g) < rank((5r)x)- 
Output: Two systems Si and 5*2 and an equation q = such that: 



The following algorithm is similar, but instead of the gcd, it returns the first input polynomial 
divided by the gcd. It is used to assimilate an inequation into a system where there already is an 
equation with the same leader, or to calculate the least common multiple of two inequations. 



i' < i 



i' < i 




where i is the quasi fiber cardinality of p and q w.r.t. S. 
Algorithm: 

1: (i,S 1 ,S 2 ) «- ResSplit (S,(S T ) X , q) 

2: (S 2 )q «- (S 2 )q U {q} 

3: return Si, S 2 , PRSi((S T )x, q, x) = 
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Algorithm 2.15 (ResSplitDivide). Input: A system S with (Sq)< x = and two polynomials p, q 
with ld(p) = ld(g) = x and a (init(p)) ^ for all a 6 &ol(S <x ). Furthermore, if rank(p) < rank(q) 
then <p B (imt(q)) ^ 0. 

Output: Two systems Si and S2 and a polynomial p such that: 

a) 5*2 = S*2 U {q} where (jSi, S 2 ^j is an i-th fibration split p and q' w.r.t. 5* 

b) p is an z-th conditional quotient of p by q' w.r.t. S 

where i is the quasi fiber cardinality of p and q' w.r.t. S, with q' = q for rank(p) > rank(<7) and 

q' = prem(<7,p, x) otherwise. 

Algorithm: 



I 


if rank(p) < rank(q) then 


2 


return ResSplitDivide(S,p, prem(<7,p, x)) 


3 


else 


4 


{i,S u S 2 ) <- ResSplit (S,p,q) 


5 


if i > then 


6 


p pquo(p, PRSi(p, prem(<7,p, x), x), x 


7 


else 


8 




9 


end if 


iO 


(S 2 ) Q ^(S 2 ) Q U{q} 


if 


return Si,S 2 ,p 


12 


end if 



Applying the last algorithm to a polynomial p and its partial derivative by its leader yields an 
algorithm to make polynomials square-free. 

In the above ResSplit-based algorithms, we had the requirement that (Sq)< x = 0. This ensures 
that all information contained in any equation of a smaller leader than x will be respected by 
reduction modulo St and thus avoids creating redundant systems. It will also be necessary for 
termination of the Decompose algorithm. This motivates the definition of a selection strategy as 
follows. 

Definition 2.16 (Select). Let Pfi n i te (M) be the set of all finite subsets of a set M. A selection 
strategy is a map 

Select : Pfi n ite({p=,P# \ p E R}) — > {p=,P^ \ p e R} : 

Q ' — >qeQ 

with the following properties: 

1. If Select(Q) = q— is an equation, then Q<i d ( 9 ) = 0- 

2. If Select(Q) = q^ is an inequation, then Q< ld ( 9 ) = 0- 

The second property of Select could be weakened further, i.e., <3<i d / g ) = 0- However, this would 
result in redundant calculations in the Decompose algorithm, thus we want all equations of the 
same leader to be treated first. 

The following algorithm is trivial. However, it will be replaced with a more complicated algo- 
rithm in §3 when the differential Thomas decomposition is treated. 

Algorithm 2.17 (InsertEquation). Input: A system S and an equation r= with ld(r) = x satisfying 
a (init(r)) ^ and cp <Xl& (r) square-free for all a G &ol(S <x ). 
Output: A system S where r= is inserted into St- 
Algorithm: 

1: if (St)x is not empty then 

2: S T ^(S T \{(S T ) X }) 

3: end if 
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4: St 4 — St U {r = } 
5: return S 

Now we present the main algorithm. It is based on all above algorithms and yields an algebraic 
Thomas decomposition. This algorithm forms the basis of the differential Thomas decomposition 
to be discussed in detail in §3. 

The general structure is as follows: In each iteration, a system S is selected from a list P of 
unfinished systems. An equation or inequation q is chosen from Sq according to the selection strat- 
egy and reduced modulo St- The algorithm assimilates q into St, potentially adding inequations 
of lower leader to Sq and adding new systems Si to P that contain a new equation of lower leader 
in (Si)q. This process works differently depending on whether q and (Sr)id(g) are equations or 
inequations, but it is based on the InitSplit, ResSplitGcd and ResSplitDivide methods in all cases. 
As soon as the algorithm yields an equation c = for c G F \ {0} or an inequation ^ in a 
system, this system is inconsistent and thus discarded. 

Algorithm 2.18 (Decompose). The algorithm is printed on page 9. 

In the next section, we consider an extension of this algorithm to partial differential systems. 
Both algorithms have been implemented, and their implementation aspects are considered in §4. 

3 The Differential Thomas Decomposition 

The differential Thomas decomposition is concerned with manipulations of polynomial differential 
equations. The basic idea for a construction of this decomposition is twofold. On the one hand a 
combinatorial calculus developed by Janet takes care of finding unique reductors and all differen- 
tial consequences by completing systems to involution. On the other hand the algebraic THOMAS 
decomposition makes the necessary splits for regularity of initials during the computation and 
ensures disjointness. 

We start by giving the basic definitions from differential algebra needed for the algorithms. 
Afterwards we summarize the Janet division and its combinatorics. The combinatorics give us a 
new algorithm InsertEquation to add equations into systems. Afterwards we review the differential 
implications of the algebraic decomposition algorithm and present the algorithm Reduce utilized 
for differential reduction. Replacing the insertion and reduction methods from the algebraic case 
with these differential methods yields the differential decomposition algorithm. 

3.1 Preliminaries from Differential Algebra 

Let A = {di, . . . ,d n } be the set of derivations (n > 0) and F be a computable Z\-differential 
field of characteristic zero. This means any dj G A is a linear operator dj : F — > F fulfilling the 
Leibniz rule. For a differential indeterminate u consider the Z\-differential polynomial ring 

F{u} := F [ Ui | i G ^>q ] , a polynomial ring infinitely generated by the algebraically independent 
set (u)a '■= {ui | i G ^>ol- The operation of dj G A on (u)a by djU\ = u\ +ej is extended linearly 
and by the Leibniz rule to F{u}. Let U = {u^\ . . . , u^} be the set of differential indeterminates. 
The multivariate Z\-differential polynomial ring is given by F{U} := F{u^} . . . {u 1 -" 1 ^}. The 

elements of {U)a '■= \ u i^ I i G ^>ci G {1,. . . , m} \ are called differential variables. 



We remark, that the algebraic closure F of F is a differential field with a differential structure 
uniquely defined by the differential structure of F (cf. [Kol73, §11.2, Lemma 1]). Let 



m 



E:=®F[[z lt ..., 



Zn}\ = F 



3 = 1 
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Algorithm 2.18 (Decompose) 

Input: A system 5" with (S') T = 0- 
Output: A Thomas decomposition of S' . 
Algorithm: 

1: P <- {S"}; Result <- 

2: while |P| > do 

3: Choose SeP;P^P\{S} 

4: if |Sq| = then 

5: Result «- EeswW U {S 1 } 

6: else 

7: g ^ Select(S Q ); S Q ^ \ {g} 

8: g <— Reduce(g, 5V); a; «— ld(g) 

9: if g £ {0^,c = | c G F \ {0}} then 

10: if x 1 then 

11: if g is an equation then 

12: if (iSt)x is an equation then 

13: if Reduce(reso((S , T)a;, q, x), St) = then 

14: (S,Si,p) «- ResSplitGCD(5,g,a;); P <- PU {Si} 

15: 5 <— lnsertEquation(S,p=) 

16: else 

17: Sq <- Sq U{q = ,res ((S T )x,q,x) = } 

18: end if 

19: else 

20: if (St)x is an inequation™ then 

21: Sq<-SqU{(St)*};St*-St\{{St)*} 

22: end if 

23: (S, S 2 ) <- lnitSplit(S, ?);?f-PU {S 2 } 

24: (S,S 3 ,p) <- ResSplitDivide(S,g, £g); P f-?U{S 3 } 

25: S <— lnsertEquation(S, p=) 

26: end if 

27: else if q is an inequation then 

28: if (St)x is an equation then 

29: (S,S 4 ,p) <- ResSplitDivide (S, (St)*, g); P <-Pu{S 4 } 

30: 5 <— lnsertEquation(S, p=) 

31: else 

32: (S, S s ) «- lnitSplit(S, g);P<-PU {S 5 } 

33: (S,S 6 ,p) <- ResSplitDivide(S,g, £g); P<-Pu{S 6 } 

34: if (St)! is an inequation then 

35: {S,S 7 ,r) <- ResSplitDivide (S, (St)*, p); P <- PU{S 7 } 

36: (St)* (r-p)* 

37: else if (St)x is empty then 

38: (St). ^ p^ 

39: end if 

40: end if 

41: end if 

42: end if 

43: P^PU{S} 

44: end if 

45: end if 

46: end while 

47: return Result 



a Remember that (St)x might be empty, and thus neither an equation nor an inequation. 
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with indeterminates z\, . . . , z n , where -F[[zi, • • • , z n ]} denotes the ring of formal power series. The 
isomorphism maps coefficients of the power series to function values of differential variables, i.e., 

a:®F[[zu...,Zn]]^ V)A :[ £ o,^,..., £ «<"> *- U h> of") 
i=l \ ieZ >o ie z >o / 

where i! := ■ . . . ■ i n \ defines the factorial of a multi-index. 

In the formulation of the algorithm the direct sum of formal power series E suffices to give a 
notion of solutions coherent to the algebraic case: For e G E we define the f-algebra homomor- 
phism 

4> e ■ F{U} -> F : ui j) h->. a(e)(u U) ) 

evaluating all differential variables of a differential polynomial at the power series e. A differential 
equation or inequation for m functions U = {u^, . . . ,u^} in n indeterminates is an element 
p G F{U} written p = or p^t, respectively. A solution of p = or p^ is an e G E with cj> e (p) = 
or 4> e (p) 7^ 0, respectively. More generally e G E is called a solution of a set P of equations 
and inequations, if it is a solution of each element in P. The set of solutions of P is denoted 
by &o\e{P) = &ol(P) C E. Since we substitute elements of F algebraically for the differential 
indeterminates, Remark (2.2), which guarantees the continuation of solutions from lower ranking 
variables to higher ranking ones, also holds here. 

Any differential F-algebra R with a differential embedding of E <— > R might be chosen as uni- 
versal set of solutions, for example a universal differential field containing F: Clearly -F[[zi, . . . , z n ]} 
embeds into its field of quotients .F((zi, . . . , z„)), and thus i^[[zi, . . . , z n ]] also embeds into a uni- 
versal differential field containing F, since F((zi, . . . , z„)) is a finitely generated differential field 
extension of F (cf. [Kol73, §11.2 and §111.7]). We denote the set of solutions in R by &ol R (P) C R. 

A finite set of equations and inequations is called a (differential) system over F{U}. We will 
be using the same notation for systems as in the algebraic Thomas decomposition introduced in 
§2.1 and §2.2, in particular a system S is represented by a pair (St, Sq). However, the candidate 
simple system St will also reflect a differential structure using combinatorial methods. We will 
elaborate on the combinatorics in the next section. 



3.2 The Combinatorics of Janet Division 

In this subsection we will focus on the combinatorics of equations, enabling us to control the infinite 
set of differential variables appearing as partial derivatives of differential indeterminates. For this 
purpose we use Janet division [GB98a] which defines these combinatorics and also automatizes 
construction of integrability conditions. An overview of modern development can be found in 
[Ger05, SeilO] and the original ideas by Janet are formulated in [Jan29]. This is achieved by 
partitioning the set of differential variables into finitely many "cones" and "free" variables. For 
creating this partition we present an algorithm for inserting new equations into an existing set 
of equations and adjusting the cone decomposition. Apart from this insertion algorithm the only 
other adaptation of the algebraic Decompose algorithm (2.18) will be the reduction algorithm 
presented in §3.3. 

We fix a (differential) ranking <, which is defined as a total order on the differential variables 
such that < dju^ and < implies djU^ < dju^ for all u^ k \u^ g U, dj G A. For 
any finite set of differential variables, a differential ranking is a ranking as defined for the algebraic 
case in §2.1. This allows us to define the largest differential variable ld(p) appearing in a differential 
polynomial p G F{U} as leader, which is set to 1 for p G F. Furthermore, define rank(p) and 
init(p) as the degree in the leader and the coefficient of \d(p) lank ( p \ respectively. Again we will 
assume 1 < for all j G {1, . . . , m} and i G Z™ . 

A set W of differential variables is closed under the action of A' C A if diiu G W V9i G A', w G 
W. The smallest such closed set containing a differential variable w denoted by (w)a' is called a 
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cone and the elements of A' we shall call (Janet) admissible derivations 4 . The Z\'-closed set 
generated by a set W of differential variables is defined to be 

(W) A > := p| Wi C (U) A . 

WiDW 

Wi Zl'-closed 

For a finite set W = {w\, . . . ,w r }, the Janet division algorithmically assigns admissible 
derivations to the elements of W such that the cones generated by the w G W are disjoint. The 
derivation di G A is assigned to the cone generated by w = G W as admissible derivation, if 
and only if 

i ; = max jij | u[P e WA' k = i k for all 1 < k < zj 

holds. We remark, that j is fixed in this definition, i.e., when constructing cones we only take 
into account other differential variables belonging to the same differential indeterminate. The 
admissible derivations assigned to w are denoted by Aw(w) C A and we call the cone (w)a w ( w ) 
the Janet cone of w with respect to W. This construction ensures disjointness of cones jDut not 
necessarily that the union of cones equals (W)a- For the Janet completion a finite set W D W 
is successively created by adding any w = diWj \S we w( w ) A w (w) to W, where Wj G W and 
di G A \ Ayy(wj). This leads to the disjoint Janet decomposition 

(W) A = tyjwu-rw 

that separates a Z\-closed set W into finitely many cones (w)a~( w ) after finitely many steps. For 
details see [Ger05, Def. 3.4] and [GB98a, Corr. 4.11]. 

With the Janet decomposition being defined for sets of differential variables, we will assign ad- 
missible derivations to differential polynomials according to their leaders. In particular, we extend 
the definitions of A\y(w) for finite W C F{U} and w G W by defining Aw(w) := Ax^f W \{\A{w)) . 

A differential polynomial q G F{U} is called reducible with respect to p G F{U}, if there exists 
i G Z^ Q such that dl 1 4 n ld(p) = ld(9{ 1 • . . . • d%p) = ld(q) and rank(aj 1 • . . . • d%p) < rank(g). 
We call a derivative of an equation by an admissible derivation an admissible prolongation. 
When restricting ourselves to admissible prolongation, we get the following concept: For a finite 
set T G F{U} we call a differential polynomial q G F{U} Janet reducible with respect to 
p G T, if there exists i G Z" such that d] 1 ■ . . . ■ d^ 1 \d(p) = ld(q) with all proper derivatives being 

admissible and rank(9 1 1 • . . . • d^p) < rank(<7) holds. We shall also say that q is Janet reducible 
modulo T if there is a p G T, such that q is Janet reducible with respect to p G T. 

A set of differential variables T C (U)a is called minimal, if for any set S C (U)a with 
l±lt€T(*)^T(t) = l±J se5 ( s )zi s ( s ) the condition T C S holds (cf. [GB98b, Def. 4.2]). We also call a set 
of differential polynomials minimal, if the corresponding set of leaders is minimal. 

In addition to the non-zero initials and square-freeness of polynomials in the candidate set St 
for a simple system (as defined in §2.2), the equations in (St) = are required to have admissible 
derivations assigned to them. When an equation p is not reducible modulo (5r) = it is added to 
(St) = and all polynomials in St with a leader being a derivative of ld(p) are removed from St, 
ensuring minimality. Furthermore, all non-admissible prolongations are created to be processed. 
This is formulated in the following algorithm: 

Algorithm 3.1 (InsertEquation). 

Input: A system 5" and a polynomial p— G F{U} not reducible modulo (S T ) = . 
Output: A system S, where (St) = C (S' t ) = U {p=} is maximal satisfying 

{ld(q) | q G (S T ) \ M} n (\d(p)) A = 0, 

S Q = S' Q U (S' T \ S T ) U {(d t q) = | q G (S T )= , di t A(St)=)(?)} • 

Algorithm: 



4 In [Ger99] and [SeilO, Chap. 7] the admissible derivations are called (Janet) multiplicative variables. 
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1: S^S' 

2: S T <r- S T U {p=} 

3: for qe S T \ {p} do 

4: if ld(q) G (ld(p)) A then 

5: Sq <- Sq U {q} 

6: S T <- S T \ {q} 

7: end if 

8: end for 

9: Reassign admissible derivations to (St) = 

10: S Q <- Sq U {(0 <g ) = | 9 G (S T ) = ,a t £ ^ (( s T )=)(«)} 

11: return S 

We remark that a non-admissible prolongations might be added to Sq again each step, even 
though it has been added before. This can be prevented by simply storing all previously generated 
non-admissible prolongations. 



3.3 Differential Simple Systems 

This section goes on reducing the differential decomposition algorithm to the algebraic one. We 
start by introducing partial solutions in order to algebraically evaluate differential polynomials at 
them yielding univariate differential polynomials. Then we present a differential reduction algo- 
rithm, as the second distinction from the algebraic decomposition algorithm. At last this section 
defines differential simple systems. 

For a differential variable x G (U)a and a power series e G E define the F- algebra homomor- 
phism 

r i — r , v i ftti^ !-> a(e)(u^), for < x 

<f> <x , e : F{U} F[v\v£ (U)a,v > x } : ^ Jj^ ^ ^ ^ 

evaluating all differential variables of a differential polynomial at e which are <-smaller than x. 

For differential reduction the Janet partition of differential variables provides the mechanism 
to get a unique reductor in a fast way (for an algorithm see [GYB01]) by restricting to admissible 
prolongations. After finding a reductor we apply a pseudo remainder algorithm (see Eq. (1)). 

We need to ensure that initials (and initials of the derivatives) of equations are non-zero. Let 
p G F{U} with x = ld(p) and define the separant sep(p) := One easily checks (cf. [Kol73, 
§1.8, lemma 5] or [Hub03b, §3.1]) that the initial of any derivative of p is sep(p) and the separant 
of any square- free equation p is nonzero on &ol(p). So by making the equations square- free, it is 
ensured that pseudo reductions are not only possible modulo p, but also modulo its derivatives. 
This allows us to formulate the differential reduction algorithm: 

Algorithm 3.2 (Reduce). 

Input: A differential system S and a polynomial p G F{U}. 

Output: A polynomial q that is not Janet reducible modulo St with 4> e {p) = if and only if 
<j> e {q) = for each e G <&ol(S). 
Algorithm: 
1: X «- ld(p) 

2: while exists q = G (St) = and ii,...,i n G Z> with ij = for dj £ A^ ST - ) =(q) such that 

dl 1 ld(g) - Id(p) and rank(^ 1 • . . . • fl*»p) > rank(q) hold do 

3: p <- prem(p, ■ . . . ■ d l n "q, x) 
4: x <— ld(p) 
5: end while 

6: if Reduce(5,init(p)) = then 

7: return Reduce(S> - init(p)a; rank (p)) 

8: else 
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9: return p 
10: end if 

A polynomial p G F{U} is called reduced 5 modulo St if Reduce^, p) = p. A polynomial 
p G F{U} reduces to q modulo St if Reduce(5 f , p) = q. 

Usually in differential algebra, one distinguishes a (full) differential reduction as used here and 
a partial (differential) reduction. Partial reduction only employs proper derivations of equations 
for reduction (cf. [Kol73, §1.9] or [Hub03b, §3.2]). This is useful for separation of differential and 
algebraic parts of the algorithm and for the use of Rosenfeld's Lemma (cf. [Ros59]). 

Definition 3.3 (Differential Simple Systems). A differential system S is (Janet) involutive, 

if all non-admissible prolongations in (St)~ reduce to zero by (Sr) = ■ 
A system S is called differentially simple or simple, if £ is 

a) algebraically simple in the finite set of differential variables appearing in it, 

b) involutive, 

c) S = is minimal, 

d) no inequation is reducible modulo S = . 

A disjoint decomposition of a system into differentially simple subsystems is called (differential) 
Thomas decomposition. 

3.4 The Differential Decomposition Algorithm 

The differential Thomas decomposition algorithm is a modification of the algebraic Thomas 
decomposition algorithm. We have already introduced the new algorithms InsertEquation (3.1) for 
adding new equations into the systems and Reduce (3.2) for reduction, that can replace their 
counterparts in the algebraic algorithm. 

Algorithm 3.4 (Differentia I Decompose). 

Input: A differential system S' with (S')t = 0- 
Output: A differential Thomas decomposition of S' . 

Algorithm: The algorithm is obtained by replacing the two subalgorithms InsertEquation and 
Reduce in (2.18) with their differential counterparts (3.1) and (3.2), respectively. 

We give an example taken from [BC99, pp. 597-600]: 

Example 3.5 (Cole-Hopf Transformation). For F := R(x,t), A = {^, ^}, and U = {rj,(} 
consider the heat equation h = % + rj xx G F{U } = and Burger's equation b = ( t + Cxx + 2(x ■ C <= 
F{U} = . To improve readability, leaders of polynomials are underlined below. 

First we claim that any power series solution for the heat equation with a non-zero constant 
term can be transformed to a solution of Burger's equation by means of the Cole-Hopf trans- 
formation A : 77 1 — >■ — . The differential Thomas decomposition for an orderly ranking with Q x > r/ t 
of 

{h = ,{r r C,-ri x ) =1 ^} 

" v ' 

consists of the single system 

S = {(vx-vO=,(vCx + Vt + vC 2 )=,!L^} 

and one checks that Reduce(S*, b) = holds. This implies that any non-zero solution of the heat 
equation is mapped by the Cole-Hopf transformation to a solution of Burger's equation. 



5 There is a fine difference between not being reducible and being reduced. In the case of not being 
reducible the initial of a polynomial can still reduce to zero and iteratively the entire polynomial. 
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In addition we claim that A is surjective. For the proof we choose an elimination ranking (cf. 
[Hub03b, §8.1] or [Bou07]) with rj 3> £, i.e., r]i > (j for all i, j G (Z>o) 2 . We compute the differential 
Thomas decomposition of {h = , b = , (r/ ■ ( — r/ x ) = , Vjt} which again consists of a single system 

S = {{Vx - V ■ C)=, (vCx+m+v C 2 )=, 6=, C_J • 

The properties of a simple system ensure that for any solution of lower ranking equations there 
exists a solution of the other equations (cf. (2.2)). The elimination ordering guarantees that the 
only constraint for £ is Burger's equation b = and thus for any solution / e &ol(b=) there exists 
a solution (g,f) € &ol(S). Furthermore, since h— was added to the input system, g e &ol(h=) 
holds and finally the equation (77 • £ — rj x ) = implies A (g) = f. 

Remark 3.6. Elements of the differential field are not subjected to splittings, unless they are 
modelled as differential indeterminates. For example to model a differential field F = C(x) with 
A = { } , we add an extra differential indeterminate X to U and replace x by X in all equations 
and inequations. We subject X to the relation J^X = 1 for X being "generic" or {-§^X — 1) • -j^X = 
0, if we allow X to degenerate to a point. This will be subject of further study. 

4 Implementation 

4.1 Algorithmic Optimizations 

In the Decompose algorithm, pseudo remainder sequences for the same pairs of polynomials are 
usually needed several times. As these calculations are expensive in general, for avoiding repeated 
calculations, it is important that the results are kept in memory and will be reused when the same 
sequence is requested again. 

If a polynomial admits factorization, we can use the it to save computation time. More precisely, 
a disjoint decomposition of the system S W {(p ■ q) = } is given by (S U {p=}, S U {p=£, q=}) and the 

system S & {(p ■ q)^} is equivalent to S U {p^ , q^ } . Let Vj := ja;j | Xj < Xi, (<S*t)^ 7^ arid := 

^Xj j Xj < Xi, (S T ) Xj = 0}. If (S T ) Xz is irreducible over the field F, := F(Zi)[Yi]/ ((S T ) <Xi ) F(z 4 )[Yi] 

for allz G {1, . . . , n}, where ((5 l T)< a;i )F(z i )[y i ] is the ideal generated by (<S l T)< a;i in the polynomial 
ring F(Zi)[Yi], factorization of polynomials can be performed over F n instead of F. 

Coefficient growth is a common problem in elimination. If possible, polynomials should be 
represented as compact as possible. Once it is known that the initial cannot vanish, the content 
(in the univariate sense) cannot vanish either. Thus, every time an initial has been added as an 
inequation to the system, one can divide the polynomial by its content. 

If the ground field F is represented as a field of fractions of a domain D (like the rationals or a 
rational function field over the rationals), it also makes sense to remove the multivariate content, 
which is an element of F. 

When reducing, in addition to reduction modulo the polynomial of the same leader, reducing 
the coefficients modulo the polynomials of lower leader can be considered. In some cases this leads 
to a reduction of sizes of coefficients, in other cases sizes increase. The latter is partly due to 
whole polynomials being multiplied with initials of the reductors. Finding a good heuristic for this 
coefficient reduction is crucial for efficiency. 

In the algebraic algorithm, polynomials don't necessarily have to be square-free when they are 
inserted into the candidate simple system. Efficiency is sometimes improved greatly by postponing 
the calculation of the square-free split as long as possible. 

In the differential case application of criteria to avoid useless reduction of non- admissible pro- 
longations can decrease computation time. The combinatorial approach used in this paper already 
avoids many reductions of so-called Z\-polynomials, as used in other approaches (see [GY06]). 
Nonetheless, using the involutive criteria 2-4 (cf. [GB98a, Ger05, AH05] and [BLOP09, §4, Prop. 
5]) which together are equivalent to the chain criterion, is valid and helpful. 

Another possible improvement is parallelization, since the main loop in line 2 of Decompose 
(2.18) can naturally be used in parallel for different systems. 
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4.2 Implementation in Maple 

Both algorithms have been implemented in the Maple computer algebra system. Packages can 
be downloaded from [BLH10], documentation and example worksheets are available there. 

The main reason for choosing Maple for the implementation is the collection of solvers for 
polynomial equations, ODEs, and PDEs already present. Furthermore, fast algorithms exist for 
polynomial factorization over finitely generated field extensions of Q and for gcd computation. 
Computation of subresultants is not available in Maple, therefore an algorithm based on [DucOO] 
is implemented for that purpose. 

Features for the differential package include arbitrary differential rankings, using functions 
implemented in Maple as differential field, computation of power series solutions, and a direct 
connection to the solvers of Maple for differential equations. 

Example 4.1. Start by loading the current version of our package: 

> with(Dif f erentialThomas) : 

> ComputeRankingC [t] , [x2,xl,y,u] , "EliminateFunction") ; 

This creates the differential polynomial ring Q{x^ 2 \ x^\ y , u} for A — {Jj}- Here u indicates 
the input, x' 1 ^ and x^ the state, and y the output of the system. The chosen ranking "<" is the 
elimination ranking with x^ 3> x^ ^>> y u, i.e., x[ 2 ^ > xj 1 ^ > j/k > ui for all i, j, k, 1 6 Z>q. 



We follow [Dio92, Ex. 1] and want to compute the external trajectories of a differential ideal 
generated by L, i.e. intersect this differential ideal with <Q{y,u}. 

> res : =Dif f erentialThomasDecomposition(L, [] ) ; 



constraints holding for lower ranking differential indeterminates can be read off the equations and 
inequations only involving these differential indeterminates, i.e., the systems shown determine the 
external trajectories of the system: 

> PrettyPrintDif f erentialSystem(res [1] ) : 

> remove (a->has (a, x2) or has (a,xl) ,%) ; 



These systems, having disjoint solution sets, are identical to the ones found in [Dio92]. 
4.3 Implementations of Similar Decomposition Algorithms 

The RegularChains package [LMX05], which is shipped with recent versions of Maple, implements 
a decomposition of a polynomial ideal into ideals represented by regular chains and a radical de- 
composition of an ideal into square-free regular chains. The solution sets of this decomposition are 
in general not disjoint. However, there is an extension called comprehensive triangular decompo- 
sition (cf. [CGL+07]) that provides disjointness on the parameters of a parametric system. The 
systems of the parameters are not simple systems though. The RegularChains package contains 
FastArithmeticTools as a subpackage implementing asymptotically fast polynomial arithmetic for 
the modular case. 



> L : = [xl [1] -u [0] *x2 [0] , x2 [1] -xl [0] -u [0] *x2 [0] , y [0] -xl [0] ] ; 

L := [.Tlx - u a x2 , x2i -xl - u x2 , y - xl ] 




[-u(t) y(t)) + (i y(t)) u(t) 2 + (i y(t)) (| u(t)) + y(t) u(t) 2 = 0, u(t) ? 0] 




16 Thomas Bachler, Vladimir Gerdt, Markus Lange-Hegermann and Daniel Robertz 

The epsilon package ([Wan03]) by Dongming Wang implements different kinds of triangular 
decompositions (including a decomposition into regular chains like above) in Maple. It is the 
only software package besides our own that implements the Thomas decomposition. It uses the 
simpler "top-down" approach that Thomas (cf. [Tho37, Tho62]) suggested, i.e., polynomials of 
higher leader are considered first. All polynomials of the same leader are combined into one common 
consequence. New systems, which contain conditions on initials of polynomials and subresultants, 
are created by splitting subalgorithms similar to ours. All these new conditions of lower leader 
are not taken into account for now and will be treated in a later step. Contrary to our approach, 
one cannot reduce modulo an unfinished system and hence inconsistency checks are less natural 
and more complicated. It is conceivable that this strategy spends too much time on computations 
with inconsistent systems. Therefore, epsilon implements highly sophisticated heuristics for early 
detection of inconsistent systems. It achieves similar performance to our implementation. 

The Maple package diffalg [BH04] deals with ordinary and partial differential equations as 
described in [BLOP09]. Its functionalities are used by symbolic differential equations solvers in 
Maple. For an input of equations and inequations it computes a radical decomposition of the 
differential ideal generated by the equations and saturated by the inequations. I.e., a description 
of the vanishing ideal of the Kolchin closure (cf. [Kol73, §IV.l]) of the solutions is computed. 
The output are differential characteristic sets as introduced by Ritt [Rit50, §1.5]. Computation 
of differential consequences is driven by reduction of Z\-polynomials, which are the analogon of 
s-polynomials in differential algebra. We found the system being optimized and well-suited for 
computations with ordinary differential equations. 

Similar algorithms as in diffalg are used in the BLAD-libraries [Bou09]. It is designed as a 
stand-alone C-library with an emphasis on usability for non-mathematicians and extensive docu- 
mentation. As it is written in C, BLAD is expected to outperform diffalg for relevant examples. 

For future publications, we plan to compare the Thomas decomposition and our implementa- 
tion with other decompositions and implementations. We also plan to further examine applications 
that benefit from the properties of simple systems. 

5 Acknowledgements 

The contents of this paper profited very much from numerous useful comments and remarks 
from Wilhclm Plesken. The authors thank him as well as Dongming Wang and Francois Boulier 
for fruitful discussions. Furthermore, our gratitude goes to the anonymous referees for valuable 
comments and for pointing out informative references. The second author (V.P.G.) acknowledges 
the Deutsche Forschungsgemeinschaft for the financial support that made his stay in Aachen 
possible. The presented results were obtained during his visits. 

References 

[AH05] Joachim Apel and Ralf Hemmecke, Detecting unnecessary reductions in an involutive ba- 
sis computation, J. Symbolic Comput. 40 (2005), no. 4-5, 1131-1149. MR MR2169107 
(2006j:13026) 14 

[BC99] Alexandru Buium and Phyllis J. Cassidy, Differential algebraic geometry and differential 
algebraic groups: from algebraic differential equations to diophantine geometry, in [Kol99] 
(1999), 567-636. 13 

[BCG+03] Y. A. Blinkov, C. F. Cid, V. P. Gerdt, W. Plesken, and D. Robertz, The MAPLE Pack- 
age Janet: I. Polynomial Systems. II. Linear Partial Differential Equations, Proc. 6th 
Int. Workshop on Computer Algebra in Scientific Computing, Passau, Germany, 2003, 
(http://wwwb.math.rwth-aachen.de/Janet), pp. 31-40 and 41-54. 2 

[BH04] Francois Boulier and Evelyne Hubert, DIFFALG: description, help pages and examples of 
use, 1996-2004, Symbolic Computation Group, University of Waterloo, Ontario, Canada 
(http : / /www- sop . inria.fr/members/Evelyne . Hubert /dif f alg/). 1, 16 

[BKRM01] Driss Bouziane, Abdelilah Kandri Rody, and Hamid Maarouf, Unmixed- dimensional decom- 
position of a finitely generated perfect differential ideal, J. Symbolic Comput. 31 (2001), no. 6, 
631-649. MR MR1834002 (2002c: 12007) 1 



Thomas Decomposition of Algebraic and Differential Systems 



17 



[BLH10] Thomas Bachler and Markus Lange-Hegermann, AlgebraicThomas and DifFeren- 
tialThomas: Thomas decomposition for algebraic and differential systems, 2008-2010, 
(http : //wwwb . math . rwth- aachen . de/thomasdecomposit ion/ ) . 15 

[BLOP95] Francois Boulier, Daniel Lazard, Francois Ollivier, and Michel Petitot, Representation for the 
radical of a finitely generated differential ideal, ISSAC, 1995, pp. 158-166. 1 

[BLOP09] , Computing representations for radicals of finitely generated differential ideals, Appl. 

Algebra Engrg. Comm. Comput. 20 (2009), no. 1, 73-121. MR MR2496662 (2010c:12005) 1, 
14, 16 

[Bou07] Francois Boulier, Differential elimination and biological modelling, Grobner bases in symbolic 
analysis, Radon Ser. Comput. Appl. Math., vol. 2, Walter de Gruyter, Berlin, 2007, pp. 109- 
137. MR MR2394771 (2009f: 12005) 14 

[Bou09] , BLAD: Bibliotheques lilloises d'algebre differentielle, 2004-2009, 

(http://www.lifl.fr/~boulier/BLAD/). 1, 16 

[CGL + 07] Changbo Chen, Oleg Golubitsky, Francois Lemaire, Marc Moreno Maza, and Wei Pan, Com- 
prehensive triangular decomposition, CASC (Victor G. Ganzha, Ernst W. Mayr, and Ev- 
genii V. Vorozhtsov, eds.), Lecture Notes in Computer Science, vol. 4770, Springer, 2007, 
pp. 73-101. 15 

[DelOO] Stephane Delliere, D.m. wang simple systems and dynamic constructible closure, Rapport de 

Recherche No. 2000-16 de l'Universite de Limoges (2000). 1 
[Dio92] Sette Diop, On universal observability, Proc. 31st Conference on Decision and Control (Tu- 

caon, Arizona), 1992. 15 

[DucOO] Lionel Ducos, Optimizations of the subresultant algorithm, J. Pure Appl. Algebra 145 (2000), 
no. 2, 149-163. MR MR1733249 (2000m:68187) 15 

[GB98a] Vladimir P. Gerdt and Yuri A. Blinkov, Involutive bases of polynomial ideals, Math. Comput. 

Simulation 45 (1998), no. 5-6, 519-541, Simplification of systems of algebraic and differential 
equations with applications. MR MR1627129 (99e:13033) 1, 2, 10, 11, 14 

[GB98b] , Minimal involutive bases, Math. Comput. Simulation 45 (1998), no. 5-6, 543- 

560, Simplification of systems of algebraic and differential equations with applications. MR 
MR1627130 (99e: 13034) 11 

[Ger99] Vladimir P. Gerdt, Completion of linear differential systems to involution, Computer alge- 
bra in scientific computing — CASC'99 (Munich), Springer, Berlin, 1999, pp. 115-137. MR 
MR1729618 (2001d:12010) 2, 10 

[Ger05] , Involutive algorithms for computing Grobner bases, Computational commutative and 

non-commutative algebraic geometry, NATO Sci. Ser. Ill Comput. Syst. Sci., vol. 196, IOS, 
Amsterdam, 2005, pp. 199-225. MR MR2179201 (2007c: 13040) 2, 10, 11, 14 

[Ger08] , On decomposition of algebraic PDE systems into simple subsystems, Acta Appl. 

Math. 101 (2008), no. 1-3, 39-51. MR MR2383543 (2009c:35003) 2 

[GY06] Vladimir P. Gerdt and Denis A. Yanovich, Investigation of the effectiveness of involutive 
criteria for computing polynomial Janet bases, Programming and Computer Software 32 
(2006), no. 3, 134-138. MR MR2267374 (2007e: 13036) 14 

[GYB01] Vladimir P. Gerdt, Denis A. Yanovich, and Yuri A. Blinkov, Fast search for the Janet divisor, 
Programming and Computer Software 27 (2001), no. 1, 22-24. MR MR1867717 12 

[Hab48] Walter Habicht, Eine Verallgemeinerung des Sturmschen Wurzelzdhlverfahrens, Comment. 
Math. Helv. 21 (1948), 99-116. MR MR0023796 (9,405f) 5 

[Hub03a] Evelyne Hubert, Notes on triangular sets and triangulation- decomposition algorithms. I. 

Polynomial systems, Symbolic and numerical scientific computation (Hagenberg, 2001), Lec- 
ture Notes in Comput. Sci., vol. 2630, Springer, Berlin, 2003, pp. 1-39. MR MR2043699 
(2005c: 13034) 1 

[Hub03b] , Notes on triangular sets and triangulation- decomposition algorithms. II. Differential 

systems, Symbolic and numerical scientific computation (Hagenberg, 2001), Lecture Notes in 
Comput. Sci., vol. 2630, Springer, Berlin, 2003, pp. 40-87. MR MR2043700 (2005c:13035) 1, 
12, 13, 14 

[Jan29] Maurice Janet, Legons sur les systemes des equationes aux derivees partielles, Cahiers Scien- 

tifiques IV, Gauthiers-Villars, Paris, 1929. 1, 10 
[Kol73] Ellis R. Kolchin, Differential algebra and algebraic groups, Academic Press, New York, 1973, 

Pure and Applied Mathematics, Vol. 54. MR MR0568864 (58 #27929) 8, 10, 12, 13, 16 
[Kol99] , Selected works of Ellis Kolchin with commentary, American Mathematical Society, 

Providence, RI, 1999, Commentaries by Armand Borel, Michael F. Singer, Bruno Poizat, 

Alexandru Buium and Phyllis J. Cassidy, Edited and with a preface by Hyman Bass, Buium 

and Cassidy. MR MR1677530 (2000g:01042) 16 



18 



Thomas Bachler, Vladimir Gerdt, Markus Lange-Hegermann and Daniel Robertz 



[LMX05] 
[LW99] 

[Mis93] 

[RiqlO] 
[Rit50] 

[Ros59] 
[SeilO] 

[sGH09] 

[Tho37] 

[Tho62] 
[Wan98] 

[WanOl] 

[Wan03] 

[Wan04] 

[WuOO] 

[YapOO] 



F. Lemaire, M. Moreno Maza, and Y. Xie, The RegularChains library in Maple, SIGSAM 
Bull. 39 (2005), no. 3, 96-97. 15 

Ziming Li and Dongming Wang, Coherent, regular and simple systems in zero decompositions 
of partial differential systems, System Science and Mathematical Sciences 12 (1999), 43-60. 
2 

Bhubaneswar Mishra, Algorithmic algebra, Texts and Monographs in Computer Science, 
Springer- Verlag, New York, 1993. MR MR1239443 (94j:68127) 5 
F. Riquier, Les systemes d'equations aux derivees partielles, 1910. 1 

Joseph F. Ritt, Differential Algebra, American Mathematical Society Colloquium Publica- 
tions, Vol. XXXIII, American Mathematical Society, New York, N. Y., 1950. MR MR0035763 
(12,7c) 1, 16 

Azriel Rosenfeld, Specializations in differential algebra, Trans. Amer. Math. Soc. 90 (1959), 
394-407. MR MR0107642 (21 #6367) 13 

Werner M. Seiler, Involution, Algorithms and Computation in Mathematics, vol. 24, Springer- 
Verlag, Berlin, 2010, The formal theory of differential equations and its applications in com- 
puter algebra. MR MR2573958 2, 10 

Xiao shan Gao and Zhenyu Huang, Efficient characteristic set algorithms for equation solv- 
ing in finite fields and application in analysis of stream ciphers, Cryptology ePrint Archive, 
Report 2009/637, 2009, http://eprint.iacr.org/. 1 

Joseph M. Thomas, Differential systems, AMS Colloquium Publications vol XXI, 1937. 1, 3, 
16 

, Systems and roots, The William Byrd Press, INC, Richmond Virginia, 1962. 1, 3, 16 

Dongming Wang, Decomposing polynomial systems into simple systems, J. Symbolic Comput. 
25 (1998), no. 3, 295-314. MR MR1615318 (99d:68130) 2, 3 

, Elimination methods, Texts and Monographs in Symbolic Computation, Springer- 

Verlag, Vienna, 2001. MR MR1826878 (2002i: 13040) 2 
, epsilon: description, help pages and 



examples of use, 



2003, 



(http://www-spiral.lip6.fr/~wang/epsilon/). 2, 16 

, Elimination practice, Imperial College Press, London, 2004, Software tools and ap- 
plications, With 1 CD-ROM (UNIX/LINUX, Windows). MR MR2050441 (2005a:68001) 2 
Wen-Tsun Wu, Mathematics mechanization, Mathematics and its Applications, vol. 489, 
Kluwer Academic Publishers Group, Dordrecht, 2000, Mechanical geometry theorem-proving, 
mechanical geometry problem-solving and polynomial equations-solving. MR MR1834540 
(2003a:01005) 1 

Chee K. Yap, Fundamental problems of algorithmic algebra, Oxford University Press, New 
York, 2000. MR MR1740761 (2000m: 12014) 5 



